A novel approach to synchronization in coupled excitable systems 
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We consider networks of coupled stochastic oscillators. When coupled we find strong collective 
oscillations, while each unit remains stochastic. In the limit N — > oo we derive a system of integro- 
delay equations and show analytically that the collective oscillations persist in a large region in 
parameter. For a regular topology with few connections between the oscillators, islands of coherent 
oscillations are formed, which merge as the amount of topological disorder increases. We link this 
transition to typical network quantities in the framework of small-world networks. 
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Collective behavior or synchronization of non- 
equilibrium systems is one of the main topics in cur- 
rent complex systems research. It can be observed in 
a variety of different physical, biological and sociologi- 
cal frameworks, such as sub-excitable media, in which 
noise- induced coherent patterns emerge 0, neurons in 
the brain, where spike patterns synchronize during ex- 
pectancy or attention |2|], or in an enthusiastic audi- 
ence which applauds in synchrony after a good perfor- 
mance ||. 

Starting from the pioneering work of Winfree Q and 
Kuramoto || ||] on coupled phase oscillators, numerous 
studies have focussed on systems where the dynamics 
of the single units is deterministic [0. Even in sys- 
tems consisting of coupled chaotic units phase locking 
may occur j8|. However, many systems are intrinsi- 
cally stochastic or driven by noise. Numerical investi- 
gations show, that synchronization can be found there as 
well H ^ 0, or that it may even be induced by 
a suitable chosen noise strength |L^, To really un- 
derstand why, and under which conditions synchroniza- 
tion appears in coupled stochastic systems, analytically 
tractable models are indispensable. There, however, lit- 
tle advancement has been achieved so far. 

In this Letter we introduce a system of N coupled dis- 
crete stochastic oscillators, each oscillator serving as a 
prototype of an excitable system. We find large coher- 
ent oscillations of the ensemble. In the limit N — > oo 
we show analytically that this coherence is due to a Hopf 
bifurcation in the dynamics of the ensemble distributions 
in case of a global coupling. Furthermore, we show nu- 
merically, that the synchronization is preserved even if 
there are only few connections between the oscillators. 
In contrast to previous investigations |T^ , [lrj| it becomes 
most pronounced for a large amount of randomness in the 
topology of the network, which we study in the frame- 
work of small world networks ||l7| . 

We consider a system of N stochastic units. Each unit 



1 * 2 

Figure 1: Three state model of an excitable unit. The process 
1 — > 2 is Markovian, while the transitions 2^3 and 3 — > 1 
are deterministic with a fixed waiting time 



is a prototypical model of a stochastic excitable oscillator 
characterized by an attracting fix point 1 from which it 
escapes by noise to an excited state 2. It then performs 
a long excursion 3 and finally returns to its rest state 1 
(Fig. p. The transition 1 — > 2 is controlled by a rate 7 
which can be expressed as an Arrhenius-like relation |fl8l| . 
The probability to stay the time t in 1 is thus given by, 



wi^2 (t) = 7cxp (—7*) 



(1) 



with mean and standard deviation I/7. The transitions 
2^3 and 3 — > 1 are deterministic with a peaked waiting 
time distributions at T2 and T3, respectively: 

w 2 ->3 (t) = S(t - r 2 ) , w 3 ^i (*) = S(t - r 3 ) . (2) 

To demonstrate that the dynamics of each individual unit 
shows typical features of excitable systems we inspect its 
power spectrum. Setting the output s(t) — 1 if the unit is 
in state 3 at time t, and s(t) — 0, otherwise, the spectrum 
is then determined employing renewal theory [Q : 
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where T = T2 + T3. The stochasticity of the process 
can be controlled by varying the amount of time spent in 
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Figure 2: Power spectrum of a single three state unit with 
t 2 = 0.1, r 3 = 1.2 and I/7 = 1.0 (solid line), I/7 = 10.0 
(dashed line). For small values of I/7 the process is oscillating 



state 1 compared to T. If I/7 3> T, the system spends 
a lot of time in state 1 and the stochasticity of the first 
step dominates the dynamics. The spectrum decreases 
then monotonically until it becomes zero at too = 2tt/t3. 
Contrary, for 1/7 <C T the process becomes oscillating. 
(Fig- 1). 

To couple the units, we first introduce the dynamic 
order parameter of the ensemble: 



/(*) 



1 N 



(4) 



with Si(t) denoting the output of unit i. A prototypical 
coupling between the individual oscillators is introduced 
by letting the rate 7 functionally depend on the value of 
f(t) in an inhibitory sigmoidal fashion: 



7(/W) = 7o 1 + Atanh 



m - r 

2(7 



(5) 



Choosing A > 0, the transition rate is large for small val- 
ues and small for large values of fit). The position and 
width of the transition is mediated by the parameters /* 
and (7. In the following the position /* is kept equal to 

0. 5. The inverse of the parameter a can be seen as an 
effective coupling parameter. If a 1, the coupling de- 
pends weakly on the value of f(t) and 7 is approximately 
equal to 70. Contrary, in the case of a <C 1, a sharp tran- 
sition between the two rates 72/1 = 70 (1 ± A) occurs if 
f(t) crosses /*. A typical trajectory of f(t) is shown in 
Fig. H for a system consisting of N — 1000 units. Al- 
though each individual unit is governed by the stochastic 
transition 1 — > 2, the whole system shows an undamped 
oscillation, which resembles coherence resonance in cou- 
pled excitable Fitz-Hugh-Nagumo units |ll], [l2^| . 

The simplicity and discrete nature of our model al- 
lows to derive an analytic condition from an integro-delay 
equation for the onset of the coherent oscillations. In the 
limit N 00 the state of the system can be described 
by the ensemble averaged occupation probabilities Pk (i) , 

1. e. the probability that a unit is in state k at time t. 
Following a mean field assumption, we identify the or- 
der parameter f(t) with Pa(t). The dynamics is then 
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Figure 3: Spike train for a globally coupled network (top) 
of N = 1000 oscillators. Each time an oscillator performs a 
2^3 transition a point is plotted. The lower panel shows 
the dynamic order parameter f(t). In all figures the following 
parameters were used: T = 2.5, T2 = 0.5, 70 = 0.5, A = 0.6, 
a = 10" 5 and /* = 0.5. 



described by the following set of integro-delay equations 



Pi(t) = 1 - P 2 (t) - P 3 (t) 

P 2 (t) = [ 7(P 3 (iO)Pi(t')^ / (6) 

m) = f T \{Pz{t'))Px{t')dt' . 

Jt-T 



While the first equation expresses the normalization con- 
dition, the second and third account for the balance of 
probability. The probabilities P2{t) and Psif) are equal 
to the time-integrated decay from state 1 from time t — T2 
up to t and from t — T up to t — T2 , respectively. Similar 
equations have been investigated before in the framework 
of delayed differential equations |^o|. For a two state rate 
equation, however, only a resonant-like behavior has been 
found ||yl . 

Now we consider the unique stationary fix-point of the 
nonlinear delay system d6h by setting Pfc(t) = Pt p2[. 
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This leads to an implicit relation for e.g. P 3 * 

P* = T3 (7\ 

3 T+1/ 7 (P 3 *)' [ '> 

which is the ratio between the time spent in state 3 and 
the mean time for one round trip. Analogous relations 
may be derived for P* and P 2 *. To investigate the sta- 
bility of the single steady state we add small perturba- 
tions P fc (t) = P£ + 5P k (t) with £fc=i<*A = 0. Using 
SPk(t) oc exp(Ai) with A ^ leads after linearization of 
Eqs. (Q) to the characteristic equation: 

1 + A" 1 | 7 (P 3 *) (1 - e' XT ) - s (e~ AT2 - e~ AT )} = 0. 

(8) 

Here we have introduced s — 7 / (P 3 *)P 1 * < and 7 '(-) as 
the first derivative of 7M. 

The solutions of Eq. (|) in A = A' + iA", A', A" € K are 
complex. A Hopf bifurcation corresponds to values where 
A crosses the imaginary axis, for which we can derive a 
condition by setting A' = 0. This gives in parametric 
dependence on the frequency A" at the bifurcation: 

7 (P 3 *) (1 - cos(A"T)) - s (cos(A"r 2 ) - cos(A"P)) = 
A" + 7 (P 3 *) sin(A"T) + s (sin(A"r 2 ) - sin(A'T)) = 0. 

Figs. ^ and |U shows the region of coherent oscillations in 
the T2 — a and 70 — a plane, respectively. For t 2 > 0, 




Figure 4: Bifurcation diagram showing the parameter region 
of coherent oscillations in the T2-U plane for different values 
of 70. The inset shows the corresponding frequency at the 
bifurcation, which is nearly independent of 70. All other pa- 
rameters as in Fig. || 

fixed 70 and a ^ there is a large region in which the en- 
semble oscillates. This region grows for increasing values 
of 70. Interestingly, there exist values of t 2 where a has 
to exceed some finite value to observe the oscillations. 
When uncoupled, the power spectra of the single units 
exhibit a maximum at finite frequency for the parame- 
ters given by the three curves shown in Fig. ||. However, 
we would like to stress that coherent oscillations of the 
ensemble can also be observed if the single uncoupled 
units are in the non-oscillating regime. The frequency at 



the bifurcation, as shown in the inset of Fig. [|, which is 
about 2.2/(27r) « 0.35 for t 2 = 0.5 agrees well with the 
frequency of the oscillations observed in the numerical 
simulation (Fig. ||) which is roughly 0.34. The depen- 
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Figure 5: Bifurcation diagram showing the parameter region 
of coherent oscillations in the 70-cr plane for different values 
of T2. All other parameters as in Fig. |i| 

dence on 70 for different values of r 2 is more complex, 
which we show in Fig. ||. Surprisingly, there are two sep- 
arated regions of coherent oscillations for small values 
of t 2 in the parameter plane (dash-dotted lines). This 
means that a path in the parameter plane in the direc- 
tion of increasing randomness, i.e. decreasing 70, may 
cause the oscillations first to vanish and then reappear. 
Upon increase of t 2 the two regions grow and merge for 
t 2 « 0.48 (solid line). The merged region maximizes at 
r 2 w 1 (dotted line) , further increase of r 2 lets it shrink 
again (dashed line). In this case, only large values of 70 
above « 14 lead to global oscillations. 

To complement the globally coupled case, we now in- 
vestigate diluted networks and the influence of the topol- 
ogy on the onset of coherent oscillations in our system. 
To study this, we consider a transition from an ordered 
topology to a random network via small-world networks, 
which have recently been introduced by Watts and Stro- 
gatz [jlTj. The transition is constructed as follows: Start- 
ing from a ring with N vertices, each vertex identified as 
one three state unit, is coupled to its k nearest neighbors 
with undirected edges. This means, that the transition 
rate of unit j now depends on the local mean field: 

m = lY, s ^> ( 9 ) 

{i,j} denoting the set of neighbors of vertex j. With 
probability p each edge is then cut and reconnected to a 
randomly chosen different vertex. In this way the param- 
eter p interpolates between a completely regular (p = 0) 
and a random network (p — 1). Please note that the 
number of connections is Nk which we require now to be 
much smaller than the number of all possible connections 
between vertices, which is N(N — l)/2. The small world 
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networks are found for values of p below « 0.1, where the 
mean shortest path between two arbitrary nodes drops 
rapidly, while the cluster index, giving the relative num- 
ber of common neighbors, is still large (see Fig. [7]). 

In Fig. [j] we show the dynamics on a completely reg- 
ular network, i.e. p — 0. There, we find small coherent 
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Figure 7: Left axis: Location of the maximum amplitude of 
the spectrum of the dynamic order parameter / (•*■). Right 
axis: Cluster index (+) and mean shortest path length (x) in 
the network. 
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Figure 6: Spike train for a regular network with p = 0, N — 
1000 and k = 18. Each time a stochastic unit undergoes a 
2^3 transition a point is plotted. Other parameters are the 
same as in Fig. [3] 

islands, which interchange over time. Since these islands 
are not in phase, there are no global oscillations. As a 
measure for the size of the global oscillations, we choose 
the height of the central peak in the power spectrum of 
f(t). For increasing randomness, i.e. p > 0, keeping k 
fixed, the oscillations become more pronounced and are 
finally maximized for complete disorder, which can be 
seen in Fig. 0. The transition to macroscopic oscillations 
shows a threshold behavior at p w 0.2. In the inset, both 
the cluster index and the mean path length are shown. 
Although there is a steep decrease in the mean path 
length for already small p 0.01, the amplitude of the 
oscillations only starts to increase, as the cluster index 
becomes smaller. This means, that in our model, com- 
pletely random networks synchronize best. For the con- 
sidered stochastic systems this stands in clear contrast 
to a previous study of chaotic systems Jl5[|, where it was 
claimed that the small-world route provides better syn- 
chronizability, compared to completely random graphs. 

In conclusion we have presented a system of TV dis- 
crete stochastic units, which constitutes a generic model 
of coupled excitable systems. When coupled we observed 
coherent oscillations. In the limit N — > 00 the system 
can be described by a system of integro-delay equations. 
A stability analysis reveals, that the fixed point of the 
dynamics becomes unstable under variation of the sys- 
tem parameters. We argue that this mechanism is valid 
for the finite system as presented by simulations. 



We further showed numerically, that the coherent os- 
cillations are preserved if the connectivity is strongly di- 
luted. Starting from a regular network with k nearest 
neighbor connections (kN <C (N — l)N/2) islands of co- 
herent oscillations with different relative phases coexist. 
With topological disorder these islands merge and macro- 
scopic oscillations of the whole ensemble can be observed. 
The transition to global oscillations is connected to the 
cluster index in the network, which drops considerably 
on the onset of the oscillations. 

The authors thank F. Wolf, M. Timme, M. Zaks and 
J. Freund for useful discussion. This work has been sup- 
ported by DFG Sfb-555. 
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